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Abstract: The energy levels of a quantum graph with time reversal sym¬ 
metry and unidirectional classical dynamics are doubly degenerate and obey 
the spectral statistics of the Gaussian Unitary Ensemble. These degenera¬ 
cies, however, are lifted when the unidirectionality is broken in one of the 
graph’s vertices by a singular perturbation. Based on a Random Matrix 
model we derive an analytic expression for the nearest neighbour distribu¬ 
tion between energy levels of such systems. As we demonstrate the result 
agrees excellently with the actual statistics for graphs with a uniform dis¬ 
tribution of eigenfunctions. Yet, it exhibits quite substantial deviations for 
classes of graphs which show strong scarring. 
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1. Introduction 

Chaotic quantum systems generically fall into three different universality classes whose 
local spectral statistics are equivalent to those of Gaussian Ensembles [l] of random 
Hermitian matrices. They are prescribed by general symmetry properties of the un¬ 
derlying Hamiltonians. In the case of a spinless, time reversal invariant ( TRI) system 
one normally expects, on the scales of mean level spacing (MLS), spectral properties 
corresponding to a Gaussian Orthogonal Ensemble ( GOE ). If TRI is broken instead a 
Gaussian Unitary Ensemble ( GUE ) describes the local spectral statistics. Finally, if the 
system has spin degrees of freedom, and TRI is not broken, its spectral statistics belong 
to the symplectic class ( GSE ). 
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Essen, Lotharstr. 1, 47048 Duisburg, Germany 
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It is well known that the above universality classification might be affected by the 
presence of additional symmetries in the system [2-4]. However, the universality class 
might also change due to special dynamical properties. In particular, this happens in 
Hamiltonians with classically unidirectional chaotic dynamics, where two directions of 
motion are ergodically separated from one another |5]-7|. Such systems can be realized 
e.g. in quantum billiards of constant width as depicted in fig. [lj The first figure (a) 
shows an experimental realisation of a billiard with smooth boundaries from |8|. In 
this case the underlying spectrum was found to consist (primarily) of quasi-degenerate 
pairs whereas they don’t obey the expected GOE statistics for TRI systems, but instead 
exhibit GUE behaviour. Qualitatively, the degeneracy can be understood from the 
fact that “left-moving” and “right-moving” modes are dynamically separated. Although 
classically these two motions are ergodically disconnected, quantum mechanically the 
two modes are still weakly coupled due to dynamical tunnelling through an integrable 
region of KAM-tori around the bouncing ball modes. 



Figure 1: (Colour online) Shown are the classical phase spaces of two constant width 
billiards with in case (a) smooth boundaries (figure taken from |8|) and in 
(b) with corners. Due to the unidirectionality of the systems trajectories in 
the upper part of the phase-space can not access the lower half, yet as the 
systems posses TRI the not shown halves look identical. On the left hand side 
the smooth boundaries give rise to bouncing ball modes (cyan) with adjacent 
KAM tori (cyan, p^O) which separate the phase space into two components, 
for further information see |8|. In contrast, the separating trajectories in the 
Reuleaux billiard, figure (b), are singular lines hitting the corners. 

For the second class, billiards with non-smooth boundaries such as the Reuleaux 
polygon shown in fig. [jjb), the splittings between quasi-degenerate states behave quite 
differently. In this case the dynamical barrier in the middle of the phase space shrinks 
to zero and the tunnelling occurs due to diffractional orbits hitting the corners of the 
billiard domain. In contrast to billiards with smooth boundaries this tunnelling effect 
is much stronger leading to large splittings comparable to the mean level spacing. See 
for instance fig. [2] showing a doublet of the Reuleaux billiard. The effect induced by the 
billiard corners is reminiscent of singular perturbations; on the semiclassical level both 
cases give rise to quite similar, singular (diffractional) classical orbits. Although it is 
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Figure 2: (Colour online) A doublet (eigenmodes 221 and 222) of the Reuleaux polygon 
whose classical phase space is shown in fig. [ljb). Dark areas correspond to 
negative values of the wavefunction, light to positive. Black “lines” in-between 
denote the nodal lines. The energy scale below is given in units of MLS. 


known, see for instance (9,10 , that normally neither the presence of corners nor singular 
perturbations affects the spectral statistics, one can anticipate essentially different results 
for unidirectional billiards. Contributions arising from diffractional orbits break down 
the unidirectionality of the billiard dynamics, thus effectively changing the universality 
class of the system. As a result, the spectral statistics of these systems do not belong to 
any of the three standard classes and it is the aim of this paper to study this phenomenon. 

Rather than considering billiards we will focus on the nearest neighbour distribution 
p u (s) of eigenvalues in quantum graphs with nearly unidirectional classical “dynam¬ 
ics”. They allow us to study the isolated effect of a singular perturbation with tunable 
strength u. Furthermore, the rank of the singular perturbation can easily be adjusted by 
prescribing the number of vertices where the unidirectionality is broken. In the present 
work we primarily study the effect of rank 1 perturbations. 


2. Unidirectional Quantum Graphs 

Quantum graphs are a widely used toy-model of chaotic quantum systems. We will first 
sketch below their general properties (for a more detailed account see for instance the 
review |TT] or 1121) and then introduce the family of unidirectional quantum graphs. 

A (closed) quantum graph consists of a set of B finite length lj edges connected at V 
vertices. The edges can be thought of as ideal ID waveguides on which a wave function 
-0 propagates. In our case the propagation is free, i.e. = — A ipj = k 2 ijjj, therefore 
'ipj = aje +lkx + bje~ lkx at each edge j. The true complexity of such systems stems from 
the boundary conditions at the graph’s vertices. At any given vertex i we are faced with 
a number of incoming waves, forming the vector , and an equal number of outgoing 
waves on the same edges forming the vector L'o,] t • The boundary condition describing 
the vertex is a unitary matrix cq matching both vectors, 

VSt = ViJl n • (!) 
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It is often referred to as “local scattering” matrix, and its unitarity ensures the conser¬ 
vation of local probability. 

Using the local scattering matrices, an internal scattering matrix S for the total graph 
can be constructed. It maps the vector of directed, incoming wave-function amplitudes 
^in = ('0fn ) ,V?n ) ) ■ ■ ■, V’in ] ) T onto the outgoing ones Cut = (vfft, m, ■ ■ •,For 
the eigenstates of the system the outgoing wave functions acquire the phases e +lklj during 
propagation along the respective edges and turn into incoming wave functions again. We 
can represent this by a self consistent equation 

S ■ diag(e ifcZl ,..., e ikl2B ) T out = T out . (2) 

V ^ > 

C(k) 


As we describe every edge by two directions, each length lj appears twice in C(k). 
The eigenenergies of the graph are now those k\ for which we can hnd a T out fulhlling 
above equation. In other words, the matrix 1 — S£(k) needs to have a zero eigenvalue 
at k — k n . If the lengths lj are non commensurate, i.e. they are given by rationally 
independent numbers, one in general expects either GOE or GUE statistics (depending 
on the symmetry of the local scatterers) for the graphs spectrum 
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Starting from this general set-up only small adjustments are needed to introduce the 
family of unidirectional quantum graphs. The classical “dynamics” on a graph T can be 
thought of as free motion of a point-like particle on edges of the graph combined with 
stochastic transitions at its vertices. The probability to pass from edge i to edge j is 
defined by the element .SVjj 2 of the scattering matrix. For a unidirectional graph we 
need to fix the vertex matrices o r in such a way that it would be impossible to switch the 
direction of motion along the edges. To this end we split an undirected graph T into two 
directed “halves” T ± such that the n-th edges of T + and T_ correspond to two possible 
directions of motions on the n-th edge of T. In addition, we require that the number of 
outgoing and incoming edges at each vertex of T± would be identical. Note that such a 
splitting of T is possible if and only if it possesses an even number of edges per vertex 
(In such a case T admits euler cycles which can be used to assign the directions along 
the edges). The following structure of the vertex scattering matrices, 


jjT o’) with Ufl} = uju, = l, (3) 

ensures that dynamics on T± are completely decoupled, see fig. [3j Here we arranged 
= ('0in + \'0in ) T such that all entries °f F + are listed first, which makes the 
underlying block structure apparent. Note that all outgoing directions on T + are incom¬ 
ing directions of T_ and the other way around. Therefore, both directions are present 
in equal number making Ui a square matrix, which is a requirement to make Ui unitary. 
Due to the off-diagonal structure of the cr’s the transition from graph T + to T_ (and 
vice versa) is impossible. In other words, a particle launched in one direction cannot 
switch to the opposite one. Further on choosing the same U t for both blocks, imply¬ 
ing Gi — aJ , ensures that the system possesses TRI. The block structure in ([3j already 
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Figure 3: (Colour online) - Right hand side: Illustration of a graph vertex with four at¬ 
tached leads which feature two separate directions, dark/blue and light/orange. 
Locally the scattering matrix does not couple those directions, e.g. waves en¬ 
tering through an incoming orange edge only leave via an outgoing orange 
edge, as indicated in Equation ([3]) . This ensures unidirectionality on a local 
scale. 

Left hand side: Directed graph T + corresponding to a unidirectional De Bruijn 
quantum graph. Its symmetric counterpart is obtained by reversing the 
direction of each edge. We highlighted two possible positions (light red/middle 
and light green/left) at which a backscattering vertex a, Equation ([5]), could 
be placed. Note though that for a rank one perturbation only one scatter can 
be introduced. We comment on the differences of both positions in Section [HJ 


indicates that the corresponding quantum graph features an exactly double degener¬ 
ate spectrum. This becomes apparent when we cast S into a block-diagonal form. To 
this end we reorder the directed edges of T such that entries of T + appear first, i.e. 
Tout = C0+\ • • •, \ • • •, Within such order S takes the form: 


S 


S 0 \ 

o s T ) ■ 


(4) 


Both directions T + and T_ form a chiral basis and are dynamically disconnected. Fur¬ 
thermore, as they are related via a transpose of S we have lj — Ij+r for the entries of the 
(also reordered) diagonal matrix C{k). 


3. Breaking Unidirectionality 


We will now consider graphs where unidirectionality is broken at just one of the graph’s 
vertices. To this end we augment one of the edges of a unidirectional graph r 0 with a 
backscattering vertex. A general 2x2 TRI backscattering matrix serving this purpose 
can be written as 


a = e ia+i7 


r 

t 



with r = ie 1 ^ sin a , t — cos a , 


(5) 


where * denotes complex conjugation. The parameter a controls the strength of back- 
scattering: In the case of a = 0 the vertex becomes transparent with the back-reflection 
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r — 0 and all incoming waves pass through (acquiring a phase e 17 ). In the opposite case, 
for a = 7r / 2 , the transmission t becomes zero and the degeneracies are lifted with the 
splittings reaching their maximum values. The parameter j3 allows back-scattered waves 
on T-t- to acquire different phases. But, the spectrum is unaffected by the choice of this 
parameter as can be inferred from semi-classical arguments. A trajectory starting on T + 
(resp. T_) on this nearly unidirectional graph needs to back-scatter twice at a to return 
to T + (T_). This leads to a cancellation of the phase associated with f3 for all closed 
orbits. The “global” phase 7 leads to a shift of the spectrum which does not change the 
spectral statistics, e.g. p„(s). 

Given such an additional scatterer the matrix S for the total graph can be written in 
the form 


S = Soil) + (e 2l “ - 1)50(7) |V>) {i>\ = £ 0 ( 7 ) ' exp ( 2i a 1^) (^D ( 6 ) 


with 


i/3/2 -i/3/2 

m = w |ta + > + ^=- |in-) . 


(7) 


Herein So is the unperturbed scattering matrix of a unidirectional quantum graph, con¬ 
taining a in its transparent form, i.e. for a — 0 , and |- 0 ) (i(j\ is a rank 1 perturbation. 
Please note that So according to our definition depends on 7 . Only for 7 = 0 the in¬ 
troduction of a does not change the spectrum of the original graph T 0 . The non-zero 
components of |in±) in eq. (j7|) correspond to the incoming directed edges on T± leading 
to the scatterer. The action of So( 7 ) 011 these vectors should be understood in the form 
So( 7 ) |in±) = e 17 |out±), i.e. as mapping incoming wave functions onto outgoing ones. 
Note that instead of adding an additional vertex to the graph we could also change the 
local scattering matrix of a 4 edge vertex of T 0 such that it corresponds to the standard 
Neumann boundary conditions (X ^ =1 ch/y = 0). In this case the decomposition (| 6 ]) of 
S into “unidirectional” and “singular” parts holds as well. In fact, all the results of 
the present paper are applicable to quantum graphs possessing this particular form of 
scattering matrix. Due to the singular type of perturbation we call these graphs nearly 
unidirectional. 

Secular Equation. We can develop some systematic insight into the graph’s spectral 
properties by looking at the eigenvalues of the unitary quantum evolution operator: 


SC(k) |A m ) = e iAm |A m ). 


( 8 ) 


Using eq. © it is straightforward to obtain a secular equation which relates the eigen¬ 
values e lAm of the perturbed system (a 7 ^ 0 ) to the doubly degenerate eigenvalues e 1£m of 
the unperturbed one (a = 0). Expanding the eigenvectors |A m ) in terms of the old a = 0 
basis |e m ) yields 


—ie _i " _ ^ e i£m \A m \ 2 
2 sin a e iA — e i£m ’ 


(9) 


where the left hand side depends only on a. Solutions to this equation in A provide 
the spectrum of eigenphases {A m } of the perturbed system. Importantly, this equation 
only determines one half of the 2 B eigenvalues of the matrix SC{k). The other half 
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are pinned to their initial values e m for a = 0. Indeed, since the original spectrum of 
the unidirectional graph is doubly degenerate one half of all eigenvalues is not affected 
by a rank 1 perturbation. (A simple way to see this is to notice that the spectra of 
the original and the perturbed graphs satisfy the interlacing property, see below.) The 
coefficients |A m | 2 = | (e m |in + ) | 2 + | (e m |in_) | 2 are the absolute values of the \e m ) eigen¬ 
vectors’ component corresponding to the edge(s) where a is located]^] In the physically 
interesting cases where one of the eigenvalues A is 0, i.e. eq. (j 2 j) holds, the corresponding 
A m are the amplitudes of the stationary wave solution on the graph. 

At first the complex secular equation appears overdetermined to give B real solutions 
for A. This, however, is not true, since the real part of this equation is trivially fulfilled 
due to the completeness condition J2m =l |A m | 2 = l, while the imaginary part yields 


B 

cota= ^2 | A m | 2 cot 

m= 1 


~2 


( 10 ) 


A prominent feature of eq. (10) is the interlacing property satisfied by its solutions, 
which is an immediate consequence of the fact that S' is a rank 1 perturbation of So- 
Specifically, for a positive a we find that the following inequality holds, e* < A, < Cj+i. 
Similarly, in the case of a < 0 we have e ,_\ < X l < This becomes apparent by noticing 
that the right hand side of equation (10), as a function of A, possesses poles at each e m . 
Therefore, the solutions A m lie in between those poles, see fig. [4} In the limiting case of 
a —» 0 they will coincide with the poles as the system’s spectrum is double degenerate. 



Figure 4: Schematic representation of the right-hand side of eqs. (10 11) as a function 
of A. One finds a solution to the secular equation whenever this sum is equal 
to (2z/)~ : , v — tan a as indicated by the horizontal, grey line. This graphical 
procedure nicely illustrates the interlacing property as all the A* can be found 
in between adjacent ej+i who define the poles of the sum. 


Large B limit. Equation (10) can be simplified further in the case of graphs with 
a large number of edges. The eigenphases e m reside on the 27T interval parametrising 
the unit circle. Their distance therefore scales like 0( 1 /b) with the number of edges. 


1 For a = 0, \A m \ 2 is identical on both edges (i.e. \ (e m |in±) | 2 = | (e m |out±) | 2 ) and independent of 

the direction of |e m ) (i.e. whether |e m ) resides on r + or T_). Yet, the incoming and outgoing 
amplitudes are related via exp (ikl^^Am^ = exp (ie m )Am Ut \ This explains the appearance of the 
phase factor in the numerator of eq. 0 - 
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Considering that cot x — 1 /x+0(x), only the N ~ \/~B closest e m will have an impact on 
the precise location of a solution A* between two e-eigenphases. To demonstrate this we 
split the sum in (10) into three parts, the central sum with N ~ \/~B elements and two 


surrounding contributions with 1 / 2 (B — N ) entries. It is straightforward to see that due 
to the periodicity of the cotangent the contributions of the outer sums on average cancel 
each other. Furthermore, it can be shown that the fluctuations of those sums scale with 
0(}/b) and therefore are negligible in the limit of large graphs. After the expansion of 
the cotangent, we obtain for the remaining sum, up to corrections of 0{}/b)\ 


J_ = ^ |4n| 2 

2 ^ i A- e m 


with v = tan a . 


( 11 ) 


The solutions of this equation obey the same interlacing property as before in eq. (10). 


Remarkably, eq. 0 is identical to a secular equation considered in [l4j. This work 
is concerned with the study of relations between old e m and new A m -eigenvalues in 
a perturbed Hamiltonian system, H = Hq + uN \ijj) (^|, where Hq is drawn from a 
GUE ensemble of N x N random matrices. Leaving aside differences in the physical 


interpretation this allows us to build upon the results of 14 in the next section, where we 


consider the spectral statistics of {A m }. As we demonstrate later, the spectral statistics 
of the graph’s physical spectrum { k n } determined by eq. (J2]) can, in the large B limit, 
be obtained from these results by a simple rescaling. 


4. Analytic Calculations. RMT approach. 

Throughout this section we derive the nearest neighbour distance distribution p u {s) for 
nearly unidirectional quantum graphs using a RMT like approach. In addition to the 
exact results presented here we provide a compact, heuristic surmise in section [ 6 } 

Recall that the perturbed spectrum consists of two parts, the pinned e-eigenvalues, 
and in between the A eigenvalues, e* < A* < ej+i, moving under perturbation (see sec. [3]). 
Because of this it is natural to split p u (s) into two parts, 

Pu(s) = \ (p?(s) + p“(s)) , (12) 

where p“(s) is the internal splitting distribution for the distances between e, and the next 
A i to the right while p® x (s) covers the expanse from A, to the next e l+ \. The derivation 
of both splitting distributions follows along the same lines but we focus mainly on p“(s) 
presenting results for p® x (s) at the end of this section. In addition, only the case v > 0 
has to be treated. As can be inferred from the secular equation, the result for v < 0 is 
obtainable by exchanging the roles of p™(s) and pjf(s), i.e. 

P?(s) = P-v(s) and p“(s) = p e *„(s). (13) 

Let us emphasize that the results we present below for v > 0 would not yield meaningful 
splitting distribution if negative v are entered, instead one should use the above relations. 







Relating the splitting distributions p(s) to the gap probability is a well established 
procedure in RMT since the later quantity is often easier to calculate, cf. (T5j. In our case 
we have two different sub-spectra and the gap probability E = E(e m ; n , e max ; A m i n , A max ) 
is dehned as probability that no eigenvalue of the respective kind, i.e. e* and A j, can be 
found in the intervals [e min , e max ] and [A min ,A max ], where, in principle, the two intervals 
should be thought of as being independent. 

To establish the connection between E and p v {s) let us first consider the probability 
to find an e-eigenvalue in an interval [e m i n — 5 e ,e m i n ], in the limit of 5 e —* 0. This is 
the probability not to have an e gap at this position and thus can be expressed as 
1 — E(e min — 5 e , e min ; 0,0). Divided by the length S e of the interval the result is the mean 
level density p ^ for the e spectrum 

P (e) (e m in) = e max ; 0,0) . (14) 

^min e m in=emax 


^min ^max 



Figure 5: (Colour online) Drawing of a gap in the {e n },{A n } spectrum with size s. To 
take derivatives of the gap probability E we expand it linearly around small 
perturbations of the respective gap boundary, which one can think of as small 
5 extensions of the original gap, see shaded intervals in the figure. 


We now extend our consideration to the gap intervals [e m i n ,e max ], [A m i n , A max ] where 
6 min = A m i n , fmax = A max while the endpoints are separated by some finite distance s = 
p^l (e min ) | A min — e max | measured on the scales of MLS p -1 , see fig. [H] Let p^C be the 
probability (density) to hnd an e-eigenvalue in the interval [e min — 5 e , e min ] simultaneously 
with another A-eigenvalue in the interval [A max , A max + hA] such that no other eigenvalues 
are present in-between. Following the same line of reasoning as before, it is given by a 
double derivative of E with respect to both edges e m i n and A max . On the other hand, 
is also equal to the probability density p ^ to hnd an e-eigenvalue in the interval 
[e min — <5 e , e min ] times the probability density to hnd the next A-eigenvalue in the distance 
s. The later one is exactly the sought p™(s). Therefore we may express it as 




d 2 


(^min ) EA m i n <9A r 


-E(e 


mm i '-max 


5 ^min? 


emin=A min =-s/2p( e )(0) ’ 

^max —Amax—H _ s/2p^ e ^ (0) 


(15) 


where, for the sake of convenience, we have chosen the gap symmetrically around zero. 
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A similar result holds for p® x (s). In principle only the boundary derivatives have to 
be exchanged, 

ex^ 1 <9 2 Jpf ^ 

Pv ( S ) = x(A) ( \ \fl\ W -^( e min, Cmax; A min , A max ) e min =A min =-s/2p( A ) (0) J V^) 

P 1 minJ £max= A max =+ s /2p( A )(0) 

where p^ A ' ) denotes the mean level density of the A-spectrum. Note that in our case 

p(A) = p( £ ). 


4.1. Gap Probability E 

We now present the derivation of the gap probability E for nearly unidirectional quan¬ 
tum graphs employing a RMT model. Our calculations rest upon the secular equation 
where the N distinct e-eigenvalues are distributed in accordance with the Circular 
Unitary Ensemble {CUE) describing the distribution of eigenvalues in TRI unidirec¬ 


tional quantum graphs. Further on, we assume, in accordance with 14 16 , a Gaussian 
distribution for the overlaps |A m | 2 : 


p(|.4m| 2 ) = N exp (-Wlylml 2 ) , 


(17) 


which can be expected for large graphs. The validity of this assumption is discussed 
towards the end of section [5] 


This RMT model was employed previously in 14 to calculate the two-point correlation 


function between the {e*} and {A*} spectra, to which end the following joint probability 
distribution P({ej}, {Aj}) for the total spectrum was derived: 


( 


p ({e*MAj}) oc 


N 


]["[ 4 sin —- - sin 


A i A j 


\ 


i,j=1 
i>j 


3 «p(-||>-«o 


(18) 


The e and A eigenvalues run from [—7r, +7 t] under the constraint of the interlacing prop¬ 
erty. Please observe that {e;} is taken without degeneracy and {A*} only consists of the 
shifted eigenvalues, similar to the usage in the secular equations. We will now utilize 


(18) in order to calculate the gap probability function: 

+7T 7r 7r +7T 

E{ 

^min? ^maxj ^min} ^max ) oc J dei JdXi Jde 2 ■ ■ ■ Jd\ N P{{c}, {Aj}) 


Ai 


€N 


N / v 

X ( ^1 Cmin)^(Cmax ^1 ^(A k A m ; n )0(A max Afc)^ J . 

k= 1 x 2 


(19) 


A special feature of the gap probability is that it relates to all eigenvalues of the spectrum 
in a uniform manner, visible in its product structure, which allows to take all integrals 
by standard methods, cf. [15]. In the following we will briefly sketch this procedure. 
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Expanding the 

sine functions into exponentials, 



N /)r 

F \ N 

N 


n sm ( e - 

-“) 0(1 II z~ l/2( - ei+ej) (e 1£i — e ie A 

= A{e i£i } JJe-s^- 1 ^, 

(20) 

i,j =1 v 

i>j 

2 ' i,o =1 

i>j 

3= 1 



one finds Vandermonde determinants A of e variables with a similar result for the A 
part. Upon reordering the domains of integration, 

+7T 7r 7r —(—7T —7T +7T +7T 62 6jV +7T 

J dei J dAi J de 2 ... dA^v ->• J dei J de 2 ... J de N J dAi ... j d\ N -i dAjv , (21) 

— 7T El Al EjV —7r El Ejv-l «1 OV-l e iV 

the integrals over the A variables can be drawn into the corresponding determinant form: 


J^dA^AO / e e 1 2 dA 1 / 2 (A 1 ) ... J^dAiMAO 

feN—l d A V— i/ 2 (AtV— l) ••• Xsjv.i dAjV— 1 /v ( A Af— 1) 

/ e ;dA^/i(A^) J^dAjv /2 (Ajv) ••• Je^dAjv/AKAjv) 

( 22 ) 

where fj(Xi) is a shorthand for exp (i(j — 1)A* + i(N — 1)A*/2 — NXi/ (2u)). The remain¬ 
ing integral over the e variables 


+7r 


+7T 


+7T 


E{e 


min? L max 


j A m ; n , A max ) oc / dei / de 2 


de 


EJV-1 


N , 

^{eJAIe^ine' 

j=i 


(23) 


can then be treated as follows. Adding the last line of the determinant D{ei} to the 
second last and then continuing recursively sets the upper boundary to 7 r in all integrals 
and shows the antisymmetry of D{ei} under exchange of e* -H- e r Therefore, the inte¬ 
grand in eq. (23) is a symmetric function, implying that the boundaries of integration 
can be extended to the full domain [—7r, +7 t] (up to a factor). 

Finally, we exploit the fact that Vandermonde determinants can be expanded into 
alternating sums 


A{e i£i } = ^(-l) H e i,j(1)ei e i,T(2)e2 


ja{N)e N 


(24) 


M 


wherein {a} represents the set of all permutations of the numbers (0, 1, ... N — 1), 
|cr| denotes the permutation’s parity and a{i) yields the number associated to i under 


permutation a. The symmetry of the integrand in (23) allows us to absorb the single 
factors into the second determinant, leading to 


-^(Cmin; Unaxj A m ; n , A max ) 


det F( ^min> £max) A m ; n , A max ) 


det F(0,0; 0, 0) 

where F(e min , e max ; A min , A max ) is a N x N matrix with elements 

+7T +7T 

F tl = j de /dAe-S( J -<) e 1 ( l -i)*-‘ i V le e 1(i - I ) i -i¥ J 


(25) 


X 


(1 - »(< 


.)«(£ 


mm ) w V max 


1 — 8(X — A min )6 ) (A max — A)^ . 


(26) 
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The denominator of E ensures that the probability to find gaps of zero width in both 
spectra is unity. It is related to the omitted normalisation of P in Equation (18). The 
integral in (26) consists of four parts, depending on the combinations of Heaviside 9- 
functions, which can be calculated explicitly. For the sake of compactness of exposition 
the resulting expressions are not presented here. They can be found in the appendix. 


4.2. Splitting Distribution 


nearest neighbour distribution as 


Using eqs. (15, 25) and taking into account that = p— N/(2n) we can write the 

(27) 


in , , _ {2 k\ d 5e d s ^etF^(sn/N,d e ,S x ) 

Pu _ 1 TV T 


S e =S X =0 


Nj d§ t det F( den )(s7r/iV, <5 e ) 
where the matrices in the numerator and denominator are given by 

F<”“>(z,Ma) := F(-z - S„z; -z,z + h) , 
F^(z,S,):=F(-z-S„-z; 0,0). 


(28) 

(29) 


For an illustration of the gap position please refer to figure [5] Expression (27) can be 
straightforwardly evaluated expanding iA num ) and iA den ) up to linear order in 5 e ■ 5 X 
(resp. S e ) and then taking the large A r -limit. (The resulting expressions are given by 
equations (53) and (54) in the appendix.) At this point it is convenient to introduce a 
pair of auxiliary TV-dimensional vectors \u) ,\u*): 

(j | u ) = e +ij ^ s , {j |u*) = e~ ij ^ s , (30) 


and N x N matrices: 


A ki= I ik + — - i 


N .N + 1 


5, 


kl 


Rki — oid — 


^sin {k — l) n /Ns ' 


(31) 


(k-l)n ) u 

Note that 1 — R is the well known sine-kernel matrix. This notation enables us to write 
the leading order expansion of both matrices (28, 29|) in a compact form: 

(32) 


F (den) = _ 2?rA + ^ + 0 ( 1 / N ) ; 


F (num) = _ 27tAR + | u ) ( u | _ g(s e A |u*) - |u*) ) ( (u\ A8 X - (u\ ) + 0(7*) , (33) 

where g = exp (i(N + l)-^s — The corrections stand here for neglected terms of 
order 0(}/n) in the elements of F. Employing the relationship 

det (A + |x) (y\) = det (A) (l + (y\ A^ 1 |x) () (34) 

the determinant of iA den ) can easily be calculated: 

det(F (den) ) = (-271-)* det(A)(l - N8 e /2ir). (35) 
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Taking also into account 


(.A + |m) {i 


-l 


= A 


-i 


A- 1 

\u) (v\ 

A" 1 

l + < 

[v\ 

A~ 

■ i 

\u) 


(36) 


we can perforin a similar expansion for the determinant of /r( num ) ; leading finally to 


p“(s) = 2 J |— det R x ( (u\ R X A | u*) — — (u\ A 1 R 1 1 u) (u\ R 1 A | u*) 
N 2 V Z7r 


+ (u\ R- 1 1 u) {u\ A-'R-'A | u*) ) + 0( 1 /n) . 

2i r J 


(37) 


Note that, although this quantity explicitly depends on N, it has a well defined limit 
for N —» oo. The leading denominator is compensated by the corresponding scaling 
of the scalar products. Furthermore, due to the presence of the exponent g the above 
expression turns out to be purely real. 


Along the same lines we are able to obtain the result for p® x (s) using eq. (16). Ex¬ 
panding the double derivative up to the leading order in 6 € ,6\ yields: 


pr« = 


det R 

~W 


((ti'l/rV) (mIa-'it'aIii) - (!i|A _1 irV) Klir'Aiti)) +o (<-/«). 

(38) 

ft is interesting to compare this result with the GUE nearest neighbour distribution 
p GUE (s) which should emerge from (38) in the limit of v —> 0. To this end note that 


p GUE (s) is related to the R matrix, 15 


, as: 


p GUE (s) = d 2 s det R = det R x (d 2 s Tr log R + {d s Tr log R ) 2 ) . 


(39) 


Further on, the derivative of R can be expressed as a rank 2 projector on |u) and \u*). 
In particular 


d s Trlogi? = Tr ((u| R 1 \u) + (u*\ R 1 \u* 


(40) 


The second derivative can be treated in the same fashion and, after some cancellation, 
we obtain 


p GUE (g) = ^(u*\R 1 \u*) (u\R 1 |w) — (u\R 1 \u*) (u*\R 


(41) 


which bears an apparent structural similarity to (38) 


5. Comparison with Quantum Graph Spectra 


Both expressions ([37J38) depend on the inverse of the sine-kernel matrix R 1 which is not 
known explicitly. Nevertheless, they serve as a very useful tool to numerically calculate 
the nearest neighbour distribution. As we neglected terms of the order 0(1/N), the 
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dimension N of the matrix R should be large enough to reach the limiting distribution 
with a sufficient precision. For practical purposes we found N = 100 to be sufficient for 
all ranges of the parameter strength u. 

The resulting nearest neighbour distributions obtained evaluating eqs. ( 37f38j ) are 
presented in fig. [6j In all cases we found a perfect agreement with the distribution 
(not shown here) of eigenvalues drawn from the Random Matrix Ensembles of sec. [4j 
Without perturbation, i.e. for v — 0, p® x (s) is given by the GUE nearest neighbour 
distribution while p™(s) is just a 5-spike at 0 due to the exact degeneracy of the system. 
Under a small perturbation this spike erodes, but both parts of the splitting distribution, 
p“(s) and p® x (s), are still distinguishable. With increasing perturbation strength both 
distributions become more similar. From the secular equation 0 one can infer that, 
in the limit v —)■ oo, they are actually identical. 

It is quite remarkable that the perturbation does not lead to a strict level repulsion 
as for any v one finds p™(0) > 0. Furthermore, under onset of the perturbation there is 
no longer a strict repulsion between \ and the e i+ \ eigenvalue to the right either, i.e. 
p® x (0) > 0 if v > 0. For comparison we present here the two point correlation function 
i? 2 (s) between the e and the A part of the spectrum which was derived in 14 m (for 
related quantities also refer to fl8|). For small distances s this function naturally agrees 
with p“(s) as is visible in fig. [bj 




s (MLS) 


Figure 6: Analytic results for p 
by 


(s), p“(s) and the Ra(s) two point correlator derived 


14 shown together for u = 2 and u = 6. As expected p“(s) and /^(s) are 
almost identical for small s. But for larger values p“(s) decays to 0 as the 
nearest neighbour can not be arbitrarily far, while /^(s) saturates towards 1 
as farther away values are no longer correlated. With increasing perturbation 
strength, see right panel, p“(s) and p® x (s) tend towards the same limiting 
distribution. 


Nearly unidirectional quantum graphs. In what follows we compare the above 
results with the nearest neighbour distributions of actual quantum graph spectra {k n }. 
For this we use several families of nearly unidirectional quantum graphs constructed 
according to the guidelines in sec. [2} To find numerical solutions to the equation 

det (l-SC(k n j) = 0 (42) 
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we used the (by definition) positive singular eigenvalues of 1— SC(k) and searched for the 
points where the lowest one becomes 0, see e.g. [19] for details of this method. Noticing 
that, for not too small graphs, the lowest lying eigenvalues depend approximately linearly 
on k it is straightforward to follow the downwards slope within a few iterations. We found 
this approach to work rather fast and accurately, such that 9 significant post-decimal 
digits can easily be achieved. 

In addition to the actual { k n } spectrum of the graph we also computed the splitting 
distribution between the eigenphases of the unitary quantum map SC(k). In the large 
B limit both spectra are known to possess the same statistics after rescaling the MLS to 
1, 11,20 . Indeed, the average “velocity” of the eigenphases over k depends on the (fixed) 
average edge length U while fluctuations decrease with increasing bond number. As a 
result, both spectral statistics coincide in the large B limit. Note, however, that since 
we obtain 2 B eigenvalues for any arbitrary k , this approach is numerically significantly 
less expensive. 

If not stated otherwise we set the back-scattering matrix d, eq. (J5|, to 


a = — arctan v , (3 = 0 , 7 = 0 

->• 


v 1 

r = - -, t = -— 


(43) 


1 — V 


1 — V 


This choice is reminiscent of placing an actual strength v h-potential on an edge ( 
H — A + uS(x)), where transmission and reflection rates of the free wave propagation 
are given by r and t. We omit here the k dependence of such type of perturbation. The 
graph lengths were chosen randomly from the interval [0,1] and rescaled such that the 
mean length is 1. 

Fully connected graph. This graph is composed of V vertices, which are interconnected 
by B — V(V — l)/2 bonds. Note that we do not allow self-loops, i.e. a vertex cannot 
be connected to itself. At each vertex the Ui scattering matrices, see eq. ([3]), are drawn 
randomly. Figs. [?] and [8] demonstrate that both p™(s) and p® x (s) agree quite nicely with 
the analytic predictions. 

Binary De Bruijn Graphs. The 2 P vertices of these graphs can be labelled by binary 
sequences of length p. Each vertex aia 2 ... a p , ai = {0,1} is connected with several 
(generically 4) others labelled by the sequences ba±a 2 ... a p _i, b = {0,1}, or a 2 a^ ... a p c, 
c = {0,1}, which are obtained by adding one symbol to the left and removing one on 
the right or vice versa, see 21 . A sketch of such a graph for p = 3 can be found in fig. [3] 
As opposed to the fully connected graphs, De Bruijn graphs feature several short cycles 
- the alternating pattern 1010 ... 0 — > 0101... 1 —> 1010 ... 0 is of length 2, while the 
vertices 11... 1 and 00... 0 have attached self-loops. The scattering matrices at each 
vertex are identical and chosen in such a way that the original graph is unidirectional. 

The results for the nearest neighbour distributions turn out to be quite sensitive to the 
choice of the back-scatterer position. If a is placed on a generic edge we once again find 
a good agreement with the RMT predictions, see fig. & ,b. As can be seen on fig. [lOj the 
same stays true if we change the scattering matrix at some generic vertex of the original 
unidirectional quantum graph T 0 to correspond to Neumann boundary conditions (thus 
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V = 1 


V = 2 




Figure 7: Histograms of the nearest neighbour distances for the 130,000 lowest k- 
eigenvalues obtained from a fully connected graph with 17 vertices plus back- 
scatterer, see text. The black line shows the analytic result derived in sec. [4] 
for comparison. The smaller dashed lines indicate the contributions of the sin¬ 
gle p™ (s) and p® x (s) to the analytical result. To distinguish them please note 
that pjf(0) > p® x (0). The lower panels demonstrate the onset of saturation 
for strong perturbations in which the distribution of both sub-splittings will 
become identical. MLS is adjusted to (e i+ i — e t ) = 1. (Weaker perturbation 
strengths can be found in fig. [8J) 




Figure 8: Histograms for p v (s) obtained from 2.7xl0 6 eigenphases of the graph’s quantum 
map. Further details as in fig. [tJ 
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breaking unidirectionality). On the other hand, if we place the back-scatterer on one of 
the self-loops (see fig. [3]) this has a drastic effect on the resulting splitting distribution, 
see fig. [9}^,d. 

To explain these findings let us recall that in the analytic RMT model we made an 
assumption, eq. (17), on the uniform (random wave) distribution of the wave-function’s 
probability density \A m \ 2 at the scatterer position. On the other hand, wave functions 
on certain types of graphs are known to exhibit enhanced localisation (scars) on some 
edges 22,23 . To shed further light onto the sensitivity of p u {s) to such enhancements 
we analysed the distribution of \A m \ 2 at different edges of the graph. Fig. 11 shows the 
results for two different edges of a De Bruijn graph illustrating the significant differences 
between the actual result and our original assumption, eq. ©■ 


“Generic” Edge “Generic” Edge 





Figure 9: Histogram for the splitting distributions in binary De Bruijn graphs. In the 
lower row the back-scattering element is located on a self-loop of the graph 
while in the upper it is far away from short cycles. It is plainly visible that the 
“self-loop” distributions deviate largely from the anticipated analytic results 
(black lines) while for the “generic” case the agreement is comparable to the 
results obtained from the fully connected graph, figs. [7] and [8j Further on, it 
appears that the effect is strongest for small splittings s and therefore mainly 
affects the internal splitting p™{s). 

Other systems. Besides the case of fixed local scattering matrices a t we also inves¬ 
tigated the case of a random choice for the De Bruijn graphs, as well as completely 
randomly constructed unidirectional graphs with low connectivity. Above deviation was 
present in all short cycles (containing up to 4 edges) of the tested graphs, decaying with 
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Figure 10: Depicted are nearest neighbour distance distributions for a De Bruijn graph 
with 64 vertices. On the left-hand side we placed, on a “generic” edge (see 
fig. 11), a general scatterer with a = 7t/4, f 3 = 7t/3 and 7 = 4/5 tt, which 
corresponds to u — 1 as indicated by the analytic black line. On the right- 
hand side no scatterer is placed at all, as a substitute we replace cq of one 
vertex with Neumann boundary conditions. The plotted analytical result 

0 the mismatch in 


represents v —> 00 . Minor deviations can be attributed 
p(\A m \ 2 ) for De Bruijn graphs, see right-hand side of fig. 
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P(U U 2 ) 


Self-Loop 



p(\A m \ 2 


“Generic” Edge 



Figure 11: Histograms showing the distribution of wave functions (absolute square) on 
two edges of a De Bruijn graph. Either on a self-loop or far away from any 
short cycles (right panel). In the generic case one finds the expected uniform 


distribution as in eq. (17). But for the self-loop the wave function either avoids 
the edge, i.e. increased probability for very small A m or it has an enhanced 
localisation. A logarithmic scaling would reveal an (exponentially decaying) 
heavy-tail for large A m where the distribution is orders of magnitudes larger 
than our assumption suggested. For instance \A m \ 2 > 0.2 denotes that at 
least one fifth of the wave-function is localised on this edge, implying a much 
smaller weight on the other 127 edges. 
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increasing cycle length. Notwithstanding, the choice of the local scattering matrices cr, 
along the cycles also has an important impact on the strength of the effect. 


6. Heuristic Surmise 


The formulas ( 37p8 ) for the nearest neighbour distribution, although exact, require to 
calculate the inverses of large matrices, ft would therefore be of interest to have a simple, 
analytical expression approximating p™(s) and p® x (s). In the following we obtain such 
an expression based on the Wigner Surmise providing the nearest neighbour distribution 
for GUE: 

" 2 / 4s 2 ^ 


p w (s) = 32 — exp- 


7r 


7r 


(44) 


First observe that independently of the perturbation strength v, the e and the shifted 
A part of the spectrum are both GUE distributed if considered separately, 10 . Yet, 


these distributions are not independent, since we have the interlacing property - a new 
eigenvalue A, is at least as far from Aj_i as from e t . To take into account correlations 
between the {e;} and {A*} spectra let us make a crude assumption that e t and A,;_i are 
separated by a fixed distance c; n . Since the distances between A*_i and A* are distributed 
according to the Wigner Surmise, the resulting distribution between e* and A*, 


p S (s , Cin) = p (s + c in )/J\f(c in ) with 


N (Qn) Q] 

7r 


+ erfc 


9 r- 


7T 


(45) 


is the sought approximation for p'^(s) (for s>0). Here A/"(c; n ) is fixed by the normali¬ 
sation condition and the optimal value of c m {y) has yet to be determined. Similarly, we 
can look at the splitting distribution p® x (s) between X t and ej+i, which we can express 
based on the same p s (s ) but with another cutting value c ex (z/). 

To identify the correct threshold c’s, we demand p s ( 0, c in ) —p'^(0), as well as p s (0, c ex ) = 
p® x (0). Using the exact solutions for s = 0 from the analytical calculations in Section |d] 
we obtain: 


K n ( 0) = 1 + - 


TrA 


-i 


p7 


( 0 ) = 1 


27 TV 


2,71V 

TrA -1 


(46) 

(47) 


see eqs. (37, 38) and, for the definition of A, (30). Although further analytic treatment 
of TrA -1 is possible, it is more convenient to use the -R 2 (s) function calculated in 14 
which is depicted in Figure [6j Recalling that i? 2 (0)=p“(0) holds we obtain 

lim Tr A -1 = 2 arctan nv . 

TV—>-0O 


(48) 


This provides us with the necessary relation to determine c in , c ex from v analytically. For 
instance, in the limit v—>oo, A/”(c in ) =p w (c in ) (likewise for c ex ) gives c in = c ex ~ 0.641. 

Surprisingly this simple surmise shows a very good agreement with the exact result 
for all ranges of perturbation strength. The comparison with the analytics is shown in 

fig- 


12 


19 
















Figure 12: Shown are the absolute deviations between the heuristic surmise p s and a 
corresponding analytical curve, which itself is plotted in the inset. In the 
cases a) and b) p™{s) is depicted for v = 0.1 and v = 1.3. The last panel 
displays p® x (s) for v = 0.3. While the error stays relatively constant for 
p“(s), in this case it decreases as v —>■ 0, whereat p® x (s) approaches the GUE 
nearest neighbour distribution. Due to the discrepancy between the exact 
nearest neighbour distribution and p w (s) the error will not vanish completely 
though. MLS is adjusted such that = 1. 


7. Conclusion and Outlook 

The main part of this paper is devoted to the analytical calculation of the nearest 
neighbour distribution p u (s) for the spectra of nearly unidirectional quantum graphs. 
Furthermore, based on the Wigner distribution for GUE we were able to obtain a sim¬ 
ple surmise giving a good approximation for p u (s) valid for an arbitrary perturbation 
strength v. These results show an excellent agreement with the data obtained from nu¬ 
merical calculations for generic ( e.g. fully connected) graphs. However, for some classes 
of graphs essential deviations were found if the perturbation is placed at edges belonging 
to short loops. It was demonstrated that such deviations can be attributed to a strong 
scarring effect at these edges. 

To investigate this scarring effect further it would be instructive to develop a semi- 
classical approach based on periodic orbit theory. This would allow the derivation of non- 
universal corrections to the RMT result based on the specific properties of the graph’s 
edges. The semi-classical approach is also needed to address the question of the splitting 
distribution in the spectrum of unidirectional billiard systems such as the Rculeaux 
polygons. As in the case of graphs, it might be expected that spectral deviations from 
standard statistics of GUE arise here clue to the presence of diffractive orbits. While 
in graphs such orbits are caused by backscattering at specially designed vertices, in the 
unidirectional billiards the same role could be played by singular classical orbits hitting 
the billiard’s corners. 

So far, our results were restricted to rank one perturbations. However, for billiards an 
effective rank of “perturbation” should (at least) depend on the number of corners. As 
the rank of the perturbation increases it is expected that the originally double degenerate 
GUE distribution must transform to the non-degenerate GOE. As a small outlook on this 
transition we numerically determined the nearest neighbour distribution for a graph with 
two/four backscattering vertices (effectively a rank two/four perturbation), presented in 
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fig. [l3j Note that, besides loosing the interlacing property and the split up into two 
sub-spectra, we find that the nearest neighbour distribution vanishes as s —> 0. This 
immediately implies level repulsion among energy levels of the system, which is absent 
for rank one perturbations. 



Figure 13: The images show the nearest neighbour distribution for a De Bruijn graph 
with 64 regular vertices and 2 (a) or 4 (b) scatterers placed on different, 
non-anomalous edges. Even for one additional scatterer (a) we see a strong 
deviation from our rank 1 result (dashed line) and instead find similarity to 
a GOE distribution (continuous line) which increases with larger numbers of 
scatterers. In both cases we find strict level repulsion as p u ( 0) = 0. 
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Appendix 


A. Calculation of Fi 


ki 


For the calculation of the double integrals in equation (26) it is convenient to expand the 


brackets and treat the four arising parts separately. Let us denote them by F^i = — 

F^ — F^ + F^ x \ where e.g. F^f stands for the integral over the 9(e — e m ; n )6 ) (e max — e) 
part only. Further on, to compactify the resulting expressions, we introduce k = k — 1 — 
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(N — l)/2 and l —l — 1 — (N — l)/2. They are given by: 
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(A) 

kl 


F, 
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sinh {njjr) - ^S~ k+I , 
111 
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fc + Z 
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6 ! 2 sillh I ~T~(A max A m ; n ) I , 

n k n l V 2 


2e n ^ 


=1 ^'/ c (emaxH-emin)/2 / _fc/ 

V 2 V 
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i^max 


2 sin ((/c + l ) 


^max ^min 


(fc+/)( £max - t“£min)/2 


nr 


fc + Z 


(49) 

(50) 


(51) 


wherein 5j. + j stands for the Kronecker-Delta which is 1 if k + l — 0 and zero otherwise 
while ni=\l — N/(2v) and n^ — ik+ N/(2v). In the case where all four ^-functions appear 
the general solution is slightly more complicated. Taking into account the position of 
the gap, see figure [5j it takes on the form: 
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kl 
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2 sin ((k + l)§s 
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k + l 
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max )/2 


n~ k ni 


e n * x min 


^ n r f m!ri '\ sinh (n 


^max ^mir 


(52) 


wherein the factor tt/N stems from the scaling of the gap boundaries. The Heaviside 
0-function distinguishes between the case of the internal splitting distribution, where 
it is one, and the external splitting distribution, where it is zero. Please observe that 


equation (52) holds only for the numerator, compare equation (27). In the case of the 


denominator F^ X> — 0. 

To obtain the expression for the numerator (in the case of internal splitting), see 


equation (27), we set e min = —sir/N — S e , A min = — sir/N , e max = +s7r/iV and A max = 
+S7r/iV + S\ A] Furthermore, we are only interested in the case where N S> 1 which 
allows for some minor simplifications. After expansion of the equations (50 52) up to 
the order 5 e 5\ we obtain: 


4° Um) - 2* (it + tt) S i+ , 


2 [il 


N\ sin ((k + l)§s 


+ 2e {il ~%W NS sinh ( ( ik + ^ s J + e 


2 v) k + l 

i(l-k+/NS-2j+ 
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~ N \ (~ N\ (~ N 


ik + + 5A 


(53) 


2 For the external splitting set: e m i n = — sir/N, A m i n = — sir/N — <5a, e m ax = stt/N + 5 e and X max = sir/N 
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In the case of the denominator we choose e m \r, = —stt/N — S e , e max = —s7r/iV and A m ; n = 
Amax = 0. Applying the same limit and expansions as for the numerator we arrive at: 


F, 


(den) 


kl 


3 i (k+V)7 


N 

27T ^ ik + 2 ^ ) S ' k+l 


+ e 


—i (I+ij'/ws 


■ r N \ r 

* + T k ■ 


(54) 


Note that the alternating sign e 1 (*’+0 7r appearing in both Fm ( 53|54 ) does not affect the 
result of the determinant and can be neglected. With the notation introduced in section 


4.2 we may cast both F^ into the matrix forms given by equations (32) and (33). 
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